To complement the fixed window analysis focused on emergence and elimination events, we also want to show that the dynamics of EWS reflect those of \(R_E\). I estimated \(R_E\) over the observations via particle filtering and state estimation. Here, I read in those values for each city and correlate them with top-performing EWS. Top-performing EWS are those with high AUC in the emergence and elimination scenarios. Note that \(R_E\) was calculated as:

\[ R_E = \frac{\beta_t}{\gamma} \frac{S_t}{N_t}. \]

library(tidyverse)
Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string
cities <- c("Agadez", "Maradi", "Niamey", "Zinder")  # the four city names
reffs <- tibble()  # empty tibble for storage
for(do_city in cities){
  fname <- paste0("../../results/filtered-states-", do_city, ".RDS")
  raw_filter <- readRDS(fname)
  reff_id <- which(raw_filter$state == "effective_r_seasonal")
  tmp_reff <- raw_filter$data[[reff_id]] %>%
    dplyr::select(time, med, week, date) %>%
    dplyr::mutate(city = do_city)
  
  reffs <- bind_rows(reffs, tmp_reff)
}
ggplot(reffs, aes(x = date, y = med, color = city)) +
  geom_hline(aes(yintercept = 1), linetype = 2, color = "grey35") +
  geom_line() +
  labs(x = "date", y = expression(italic(R)[E]), 
       title = "Effective reproduction number over time") +
  ggthemes::scale_color_colorblind()

Now I calculate the EWS at each observation time using the spaero::get_stats() function. I use a bandwidth of 30 weeks and set backward_only = TRUE so that only past data are used when calculating EWS.

library(spaero)
fname <- "../../data/clean-data/weekly-measles-incidence-niger-cities-clean.RDS"
measles_data <- readRDS(fname) %>%
  filter(year > 1994)  # drop first NA year, only used for modeling
all_stats <- tibble()  # empty tibble for storage
for(do_region in unique(measles_data$region)){
  
  cases <- measles_data %>%
    filter(region == do_region) %>%
    pull(cases)
  
  city_stats <- spaero::get_stats(
    x = cases,
    center_trend = "local_constant", 
    center_kernel = "uniform", 
    center_bandwidth = 30, 
    stat_trend = "local_constant", 
    stat_kernel = "uniform", 
    stat_bandwidth = 30, 
    lag = 1, 
    backward_only = TRUE
  )$stats
  
  city_stats_tb <- as_tibble(city_stats) %>%
    mutate(
      time_iter = 1:n(),
      date = unique(measles_data$date)
    ) %>%
    gather(key = ews, value = value, -time_iter, -date) %>%
    mutate(region = do_region)
  
  all_stats <- bind_rows(all_stats, city_stats_tb)
}
# best_ews <- c("autocorrelation", "autocovariance", "mean", "variance")
# all_stats <- all_stats %>%
#   dplyr::filter(ews %in% best_ews)

Now I look at scatterplots of \(R_E\) versus EWS.

reffs <- reffs %>%
  dplyr::mutate(region = paste0(city, " (City)")) %>%
  dplyr::select(-city)
ews_reffs <- left_join(all_stats, reffs)
Joining, by = c("date", "region")
ews_reffs_subcrit <- ews_reffs %>%
  filter(med < 1)
ggplot(ews_reffs_subcrit, aes(x = value, y = med, color = region)) +
  geom_point(alpha = 0.2) +
  facet_wrap(~ews, scales = "free") +
  labs(x = "EWS value", y = expression(italic(R)[E])) +
  ggthemes::scale_color_colorblind()

# for(do_region in unique(ews_reffs$region)){
#   print(ggplot(filter(ews_reffs, region == do_region & ews == "variance"), 
#        aes(x = med, y = value, color = lubridate::week(date))) +
#     geom_point() +
#     facet_wrap(~lubridate::year(date), scales = "free") +
#     labs(x = expression(R[E]), y = "EWS value", title = do_region) +
#     scale_color_viridis_c(name = "week of year", direction = 1))
# }

No clear signal in those plots. Next I’ll look at the first difference of each EWS (x) between time t and \(t+1\). I calculate the difference in three ways:

  1. Simple difference: \(x(t) - x(t-1)\)
  2. Ratio: \(\frac{x(t)}{x(t-1)}\)
  3. Log ratio: \(\text{log}\left(\frac{x(t)}{x(t-1)}\right)\)
ews_reffs_diffed <- ews_reffs %>%
  group_by(ews, region) %>%
  mutate(
    diff_value = value - lag(value, 1),
    ratio_value = value / lag(value, 1),
    log_ratio_value = log(ratio_value)
  ) %>%
  rename("Reff" = med)
NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

Correlations between \(R_E\) and first differenced EWS

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  ggplot(aes(y = Reff, x = diff_value)) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(x = "First difference of EWS", y = expression(R[E])) +
  theme(legend.position = "top")

Correlations between \(R_E\) and “first ratioed” EWS

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  ggplot(aes(y = Reff, x = ratio_value)) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(x = "First ratio of EWS", y = expression(R[E])) +
  theme(legend.position = "top")

Correlations between \(R_E\) and log of “first ratioed” EWS

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  ggplot(aes(y = Reff, x = log_ratio_value)) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(x = "log(First ratio of EWS)", y = expression(R[E])) +
  theme(legend.position = "top")

Correlations between \(R_E\) and log of “first ratioed” EWS, near-zeros removed

It looks like there might be trends when the EWS are actually “moving”, that is, when \(\text{log}\left(\frac{x(t)}{x(t-1)} \right) \ne 0\) or approximately so. So, I’ll remove those values that are close to zero by assuming all values less than \(|0.05|\) are essentially zero. Now some patterns emerge as expected.

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  filter(round(log_ratio_value,1) != 0) %>%
  ggplot(aes(x = log_ratio_value, y = Reff)) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  stat_smooth(method = "lm", se = FALSE, linetype = 1, color = "black") +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(x = "log(First ratio of EWS)", y = expression(R[E])) +
  theme(legend.position = "top")

Calculate correlations after removing ~0 EWS diffs

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  filter(round(log_ratio_value,1) != 0) %>%
  group_by(ews, region) %>%
  summarise(corr_with_reff = cor(log_ratio_value, Reff, 
                                 use = "pairwise.complete.obs", 
                                 method = "spearman")) %>%
  ggplot(aes(x = ews, y = corr_with_reff, fill = region)) +
  geom_col(position = position_dodge(), width = 0.5) +
  geom_hline(aes(yintercept = 0), col = "grey40") +
  labs(y = "Spearman's rank correlation", x = "EWS") +
  ggthemes::scale_fill_colorblind() +
  scale_y_continuous(limits = c(-0.6,0.6)) +
  coord_flip()

testdf <- ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  filter(round(log_ratio_value,1) != 0) %>%
  filter(ews == "variance")
summary(lm(Reff~log_ratio_value, data = testdf))

Call:
lm(formula = Reff ~ log_ratio_value, data = testdf)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46769 -0.12997 -0.00631  0.16118  0.45138 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.746550   0.008797  84.866  < 2e-16 ***
log_ratio_value 0.248097   0.030205   8.214 1.77e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1753 on 510 degrees of freedom
Multiple R-squared:  0.1168,    Adjusted R-squared:  0.1151 
F-statistic: 67.47 on 1 and 510 DF,  p-value: 1.775e-15

Look at lag in dynamics

Due to the stochastic seasonal dynamics, there may be a lag between \(R_E\) and EWS that are calculated from raw case data. To look at this, I will plot \(R_E\) and each EWS on the same plot, with peaks in each year highlighted. These plots will come from simulated series where we know that the \(R_E\) and the dynamics are explicitly linked (unlike \(R_E\) estimates from particle filtering with the real data).

all_sim_files <- list.files("../../simulations/")
sim_ids <- grep("bootstrap*", all_sim_files)
sim_files <- all_sim_files[sim_ids]
all_sims <- tibble()  # empty tibble for storage
for(do_file in sim_files){
  tmp_sim <- readRDS(paste0("../../simulations/", do_file)) %>%
    filter(sim == 1) %>%
    unnest() %>%
    mutate(city = str_sub(do_file, 16, -5))
  
  all_sims <- bind_rows(all_sims, tmp_sim)
}
sim_ews <- tibble()  # empty storage df
for(do_city in unique(all_sims$city)){
  tmp_data <- all_sims %>%
    filter(city == do_city)
  
  tmp_cases <- tmp_data$reports
  
  tmp_stats <- spaero::get_stats(
    x = tmp_cases,
    center_trend = "local_constant", 
    center_kernel = "uniform", 
    center_bandwidth = 30, 
    stat_trend = "local_constant", 
    stat_kernel = "uniform", 
    stat_bandwidth = 30, 
    lag = 1, 
    backward_only = TRUE
  )$stats
  
  sim_ews <- sim_ews %>%
    bind_rows(
      tmp_data %>% bind_cols(as_tibble(tmp_stats))
    )
}
sim_ews <- sim_ews %>%
  mutate(critical = ifelse(RE_seas < 1, FALSE, TRUE))
for(do_city in unique(sim_ews$city)){
  tmp <- filter(sim_ews, city == do_city)
  x <- tmp$time
  y <- tmp$RE_seas
  y2 <- log(tmp$variance)
  color <- ifelse(tmp$critical == TRUE, "red", "blue")
  
  par(mar = c(5,5,2,5))
  plot(x, y, type = "n", ylab = expression(italic(R)[E]), xlab = "Date", main = do_city)
  segments(x0 = x, y0 = y, x1 = x+diff(x), y1 = y+diff(y), col = color)
  par(new = T)
  plot(x, y2, type = "l", axes=F, xlab=NA, ylab=NA, lwd = 1.5)
  axis(side = 4)
  mtext(side = 4, line = 3, "Log variance")
}

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

for(do_city in unique(sim_ews$city)){
  tmp <- filter(sim_ews, city == do_city)
  x <- tmp$time
  y <- tmp$RE_seas
  y2 <- tmp$autocorrelation
  color <- ifelse(tmp$critical == TRUE, "red", "blue")
  
  par(mar = c(5,5,2,5))
  plot(x, y, type = "n", ylab = expression(italic(R)[E]), xlab = "Date", main = do_city)
  segments(x0 = x, y0 = y, x1 = x+diff(x), y1 = y+diff(y), col = color)
  par(new = T)
  plot(x, y2, type = "l", axes=F, xlab=NA, ylab=NA, lwd = 1.5)
  axis(side = 4)
  mtext(side = 4, line = 3, "Autocorrelation")
}

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

for(do_city in unique(sim_ews$city)){
  tmp <- filter(sim_ews, city == do_city)
  x <- tmp$time
  y <- tmp$RE_seas
  y2 <- log(tmp$decay_time + 1)
  color <- ifelse(tmp$critical == TRUE, "red", "blue")
  
  par(mar = c(5,5,2,5))
  plot(x, y, type = "n", ylab = expression(italic(R)[E]), xlab = "Date", main = do_city)
  segments(x0 = x, y0 = y, x1 = x+diff(x), y1 = y+diff(y), col = color)
  par(new = T)
  plot(x, y2, type = "l", axes=F, xlab=NA, ylab=NA, lwd = 1.5)
  axis(side = 4)
  mtext(side = 4, line = 3, "Log decay time (+1)")
}

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

LS0tCnRpdGxlOiAiQ29ycmVsYXRpb25zIGJldHdlZW4gRVdTIGFuZCBlZmZlY3RpdmUgcmVwcm9kdWN0aW9uIG51bWJlciIKb3V0cHV0OiBodG1sX25vdGVib29rCmRhdGU6IDIwMTktMDEtMjIKLS0tCgpUbyBjb21wbGVtZW50IHRoZSBmaXhlZCB3aW5kb3cgYW5hbHlzaXMgZm9jdXNlZCBvbiBlbWVyZ2VuY2UgYW5kIGVsaW1pbmF0aW9uIGV2ZW50cywgd2UgYWxzbyB3YW50IHRvIHNob3cgdGhhdCB0aGUgZHluYW1pY3Mgb2YgRVdTIHJlZmxlY3QgdGhvc2Ugb2YgJFJfRSQuCkkgZXN0aW1hdGVkICRSX0UkIG92ZXIgdGhlIG9ic2VydmF0aW9ucyB2aWEgcGFydGljbGUgZmlsdGVyaW5nIGFuZCBzdGF0ZSBlc3RpbWF0aW9uLgpIZXJlLCBJIHJlYWQgaW4gdGhvc2UgdmFsdWVzIGZvciBlYWNoIGNpdHkgYW5kIGNvcnJlbGF0ZSB0aGVtIHdpdGggdG9wLXBlcmZvcm1pbmcgRVdTLgpUb3AtcGVyZm9ybWluZyBFV1MgYXJlIHRob3NlIHdpdGggaGlnaCBBVUMgaW4gdGhlIGVtZXJnZW5jZSBhbmQgZWxpbWluYXRpb24gc2NlbmFyaW9zLgpOb3RlIHRoYXQgJFJfRSQgd2FzIGNhbGN1bGF0ZWQgYXM6CgokJApSX0UgPSBcZnJhY3tcYmV0YV90fXtcZ2FtbWF9IFxmcmFje1NfdH17Tl90fS4KJCQKCmBgYHtyIGxvYWQtcmV9CmxpYnJhcnkodGlkeXZlcnNlKQoKY2l0aWVzIDwtIGMoIkFnYWRleiIsICJNYXJhZGkiLCAiTmlhbWV5IiwgIlppbmRlciIpICAjIHRoZSBmb3VyIGNpdHkgbmFtZXMKcmVmZnMgPC0gdGliYmxlKCkgICMgZW1wdHkgdGliYmxlIGZvciBzdG9yYWdlCgpmb3IoZG9fY2l0eSBpbiBjaXRpZXMpewogIGZuYW1lIDwtIHBhc3RlMCgiLi4vLi4vcmVzdWx0cy9maWx0ZXJlZC1zdGF0ZXMtIiwgZG9fY2l0eSwgIi5SRFMiKQogIHJhd19maWx0ZXIgPC0gcmVhZFJEUyhmbmFtZSkKICByZWZmX2lkIDwtIHdoaWNoKHJhd19maWx0ZXIkc3RhdGUgPT0gImVmZmVjdGl2ZV9yX3NlYXNvbmFsIikKICB0bXBfcmVmZiA8LSByYXdfZmlsdGVyJGRhdGFbW3JlZmZfaWRdXSAlPiUKICAgIGRwbHlyOjpzZWxlY3QodGltZSwgbWVkLCB3ZWVrLCBkYXRlKSAlPiUKICAgIGRwbHlyOjptdXRhdGUoY2l0eSA9IGRvX2NpdHkpCiAgCiAgcmVmZnMgPC0gYmluZF9yb3dzKHJlZmZzLCB0bXBfcmVmZikKfQoKZ2dwbG90KHJlZmZzLCBhZXMoeCA9IGRhdGUsIHkgPSBtZWQsIGNvbG9yID0gY2l0eSkpICsKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gMSksIGxpbmV0eXBlID0gMiwgY29sb3IgPSAiZ3JleTM1IikgKwogIGdlb21fbGluZSgpICsKICBsYWJzKHggPSAiZGF0ZSIsIHkgPSBleHByZXNzaW9uKGl0YWxpYyhSKVtFXSksIAogICAgICAgdGl0bGUgPSAiRWZmZWN0aXZlIHJlcHJvZHVjdGlvbiBudW1iZXIgb3ZlciB0aW1lIikgKwogIGdndGhlbWVzOjpzY2FsZV9jb2xvcl9jb2xvcmJsaW5kKCkKYGBgCgpOb3cgSSBjYWxjdWxhdGUgdGhlIEVXUyBhdCBlYWNoIG9ic2VydmF0aW9uIHRpbWUgdXNpbmcgdGhlIGBzcGFlcm86OmdldF9zdGF0cygpYCBmdW5jdGlvbi4KSSB1c2UgYSBiYW5kd2lkdGggb2YgMzAgd2Vla3MgYW5kIHNldCBgYmFja3dhcmRfb25seSA9IFRSVUVgIHNvIHRoYXQgb25seSBwYXN0IGRhdGEgYXJlIHVzZWQgd2hlbiBjYWxjdWxhdGluZyBFV1MuCgpgYGB7ciBld3MtY2FsY30KbGlicmFyeShzcGFlcm8pCgpmbmFtZSA8LSAiLi4vLi4vZGF0YS9jbGVhbi1kYXRhL3dlZWtseS1tZWFzbGVzLWluY2lkZW5jZS1uaWdlci1jaXRpZXMtY2xlYW4uUkRTIgptZWFzbGVzX2RhdGEgPC0gcmVhZFJEUyhmbmFtZSkgJT4lCiAgZmlsdGVyKHllYXIgPiAxOTk0KSAgIyBkcm9wIGZpcnN0IE5BIHllYXIsIG9ubHkgdXNlZCBmb3IgbW9kZWxpbmcKCmFsbF9zdGF0cyA8LSB0aWJibGUoKSAgIyBlbXB0eSB0aWJibGUgZm9yIHN0b3JhZ2UKCmZvcihkb19yZWdpb24gaW4gdW5pcXVlKG1lYXNsZXNfZGF0YSRyZWdpb24pKXsKICAKICBjYXNlcyA8LSBtZWFzbGVzX2RhdGEgJT4lCiAgICBmaWx0ZXIocmVnaW9uID09IGRvX3JlZ2lvbikgJT4lCiAgICBwdWxsKGNhc2VzKQogIAogIGNpdHlfc3RhdHMgPC0gc3BhZXJvOjpnZXRfc3RhdHMoCiAgICB4ID0gY2FzZXMsCiAgICBjZW50ZXJfdHJlbmQgPSAibG9jYWxfY29uc3RhbnQiLCAKICAgIGNlbnRlcl9rZXJuZWwgPSAidW5pZm9ybSIsIAogICAgY2VudGVyX2JhbmR3aWR0aCA9IDMwLCAKICAgIHN0YXRfdHJlbmQgPSAibG9jYWxfY29uc3RhbnQiLCAKICAgIHN0YXRfa2VybmVsID0gInVuaWZvcm0iLCAKICAgIHN0YXRfYmFuZHdpZHRoID0gMzAsIAogICAgbGFnID0gMSwgCiAgICBiYWNrd2FyZF9vbmx5ID0gVFJVRQogICkkc3RhdHMKICAKICBjaXR5X3N0YXRzX3RiIDwtIGFzX3RpYmJsZShjaXR5X3N0YXRzKSAlPiUKICAgIG11dGF0ZSgKICAgICAgdGltZV9pdGVyID0gMTpuKCksCiAgICAgIGRhdGUgPSB1bmlxdWUobWVhc2xlc19kYXRhJGRhdGUpCiAgICApICU+JQogICAgZ2F0aGVyKGtleSA9IGV3cywgdmFsdWUgPSB2YWx1ZSwgLXRpbWVfaXRlciwgLWRhdGUpICU+JQogICAgbXV0YXRlKHJlZ2lvbiA9IGRvX3JlZ2lvbikKICAKICBhbGxfc3RhdHMgPC0gYmluZF9yb3dzKGFsbF9zdGF0cywgY2l0eV9zdGF0c190YikKfQoKIyBiZXN0X2V3cyA8LSBjKCJhdXRvY29ycmVsYXRpb24iLCAiYXV0b2NvdmFyaWFuY2UiLCAibWVhbiIsICJ2YXJpYW5jZSIpCiMgYWxsX3N0YXRzIDwtIGFsbF9zdGF0cyAlPiUKIyAgIGRwbHlyOjpmaWx0ZXIoZXdzICVpbiUgYmVzdF9ld3MpCmBgYAoKTm93IEkgbG9vayBhdCBzY2F0dGVycGxvdHMgb2YgJFJfRSQgdmVyc3VzIEVXUy4KCmBgYHtyIHNjYXR0ZXJzfQpyZWZmcyA8LSByZWZmcyAlPiUKICBkcGx5cjo6bXV0YXRlKHJlZ2lvbiA9IHBhc3RlMChjaXR5LCAiIChDaXR5KSIpKSAlPiUKICBkcGx5cjo6c2VsZWN0KC1jaXR5KQoKZXdzX3JlZmZzIDwtIGxlZnRfam9pbihhbGxfc3RhdHMsIHJlZmZzKQpld3NfcmVmZnNfc3ViY3JpdCA8LSBld3NfcmVmZnMgJT4lCiAgZmlsdGVyKG1lZCA8IDEpCgpnZ3Bsb3QoZXdzX3JlZmZzX3N1YmNyaXQsIGFlcyh4ID0gdmFsdWUsIHkgPSBtZWQsIGNvbG9yID0gcmVnaW9uKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICBmYWNldF93cmFwKH5ld3MsIHNjYWxlcyA9ICJmcmVlIikgKwogIGxhYnMoeCA9ICJFV1MgdmFsdWUiLCB5ID0gZXhwcmVzc2lvbihpdGFsaWMoUilbRV0pKSArCiAgZ2d0aGVtZXM6OnNjYWxlX2NvbG9yX2NvbG9yYmxpbmQoKQoKIyBmb3IoZG9fcmVnaW9uIGluIHVuaXF1ZShld3NfcmVmZnMkcmVnaW9uKSl7CiMgICBwcmludChnZ3Bsb3QoZmlsdGVyKGV3c19yZWZmcywgcmVnaW9uID09IGRvX3JlZ2lvbiAmIGV3cyA9PSAidmFyaWFuY2UiKSwgCiMgICAgICAgIGFlcyh4ID0gbWVkLCB5ID0gdmFsdWUsIGNvbG9yID0gbHVicmlkYXRlOjp3ZWVrKGRhdGUpKSkgKwojICAgICBnZW9tX3BvaW50KCkgKwojICAgICBmYWNldF93cmFwKH5sdWJyaWRhdGU6OnllYXIoZGF0ZSksIHNjYWxlcyA9ICJmcmVlIikgKwojICAgICBsYWJzKHggPSBleHByZXNzaW9uKFJbRV0pLCB5ID0gIkVXUyB2YWx1ZSIsIHRpdGxlID0gZG9fcmVnaW9uKSArCiMgICAgIHNjYWxlX2NvbG9yX3ZpcmlkaXNfYyhuYW1lID0gIndlZWsgb2YgeWVhciIsIGRpcmVjdGlvbiA9IDEpKQojIH0KYGBgCgpObyBjbGVhciBzaWduYWwgaW4gdGhvc2UgcGxvdHMuCk5leHQgSSdsbCBsb29rIGF0IHRoZSBmaXJzdCBkaWZmZXJlbmNlIG9mIGVhY2ggRVdTICgqeCopIGJldHdlZW4gdGltZSAqdCogYW5kICR0KzEkLgpJIGNhbGN1bGF0ZSB0aGUgZGlmZmVyZW5jZSBpbiB0aHJlZSB3YXlzOgoKMS4gU2ltcGxlIGRpZmZlcmVuY2U6ICR4KHQpIC0geCh0LTEpJAoyLiBSYXRpbzogJFxmcmFje3godCl9e3godC0xKX0kCjMuIExvZyByYXRpbzogJFx0ZXh0e2xvZ31cbGVmdChcZnJhY3t4KHQpfXt4KHQtMSl9XHJpZ2h0KSQKCmBgYHtyIGNhbGMtZXdzLWRpZmZzfQpld3NfcmVmZnNfZGlmZmVkIDwtIGV3c19yZWZmcyAlPiUKICBncm91cF9ieShld3MsIHJlZ2lvbikgJT4lCiAgbXV0YXRlKAogICAgZGlmZl92YWx1ZSA9IHZhbHVlIC0gbGFnKHZhbHVlLCAxKSwKICAgIHJhdGlvX3ZhbHVlID0gdmFsdWUgLyBsYWcodmFsdWUsIDEpLAogICAgbG9nX3JhdGlvX3ZhbHVlID0gbG9nKHJhdGlvX3ZhbHVlKQogICkgJT4lCiAgcmVuYW1lKCJSZWZmIiA9IG1lZCkKYGBgCgojIyMgQ29ycmVsYXRpb25zIGJldHdlZW4gJFJfRSQgYW5kIGZpcnN0IGRpZmZlcmVuY2VkIEVXUwoKYGBge3IgZmlyc3QtZGlmZi1zY2F0dGVycywgZmlnLmhlaWdodD0yMCwgZmlnLndpZHRoPTh9CmV3c19yZWZmc19kaWZmZWQgJT4lCiAgZmlsdGVyKFJlZmYgPCAxKSAlPiUKICBnZ3Bsb3QoYWVzKHkgPSBSZWZmLCB4ID0gZGlmZl92YWx1ZSkpICsKICBnZW9tX3BvaW50KGFlcyhjb2xvciA9IGx1YnJpZGF0ZTo6d2VlayhkYXRlKSksIGFscGhhID0gMC43KSArCiAgZmFjZXRfd3JhcChld3N+cmVnaW9uLCBzY2FsZXMgPSAiZnJlZSIsIG5jb2wgPSA0KSArCiAgc2NhbGVfY29sb3JfdmlyaWRpc19jKG5hbWUgPSAiV2VlayBvZiB5ZWFyIiwgbGltaXRzID0gYygxLDUyKSwgCiAgICAgICAgICAgICAgICAgICAgICAgIGJyZWFrcyA9IGMoMSwyNiw1MikpICsKICBsYWJzKHggPSAiRmlyc3QgZGlmZmVyZW5jZSBvZiBFV1MiLCB5ID0gZXhwcmVzc2lvbihSW0VdKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiKQpgYGAKCiMjIyBDb3JyZWxhdGlvbnMgYmV0d2VlbiAkUl9FJCBhbmQgImZpcnN0IHJhdGlvZWQiIEVXUwoKYGBge3IgZmlnLmhlaWdodD0yMCwgZmlnLndpZHRoPTh9CmV3c19yZWZmc19kaWZmZWQgJT4lCiAgZmlsdGVyKFJlZmYgPCAxKSAlPiUKICBnZ3Bsb3QoYWVzKHkgPSBSZWZmLCB4ID0gcmF0aW9fdmFsdWUpKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3IgPSBsdWJyaWRhdGU6OndlZWsoZGF0ZSkpLCBhbHBoYSA9IDAuNykgKwogIGZhY2V0X3dyYXAoZXdzfnJlZ2lvbiwgc2NhbGVzID0gImZyZWUiLCBuY29sID0gNCkgKwogIHNjYWxlX2NvbG9yX3ZpcmlkaXNfYyhuYW1lID0gIldlZWsgb2YgeWVhciIsIGxpbWl0cyA9IGMoMSw1MiksIAogICAgICAgICAgICAgICAgICAgICAgICBicmVha3MgPSBjKDEsMjYsNTIpKSArCiAgbGFicyh4ID0gIkZpcnN0IHJhdGlvIG9mIEVXUyIsIHkgPSBleHByZXNzaW9uKFJbRV0pKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpCmBgYAoKIyMjIENvcnJlbGF0aW9ucyBiZXR3ZWVuICRSX0UkIGFuZCBsb2cgb2YgImZpcnN0IHJhdGlvZWQiIEVXUwoKYGBge3IgZmlnLmhlaWdodD0yMCwgZmlnLndpZHRoPTh9CmV3c19yZWZmc19kaWZmZWQgJT4lCiAgZmlsdGVyKFJlZmYgPCAxKSAlPiUKICBnZ3Bsb3QoYWVzKHkgPSBSZWZmLCB4ID0gbG9nX3JhdGlvX3ZhbHVlKSkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG9yID0gbHVicmlkYXRlOjp3ZWVrKGRhdGUpKSwgYWxwaGEgPSAwLjcpICsKICBmYWNldF93cmFwKGV3c35yZWdpb24sIHNjYWxlcyA9ICJmcmVlIiwgbmNvbCA9IDQpICsKICBzY2FsZV9jb2xvcl92aXJpZGlzX2MobmFtZSA9ICJXZWVrIG9mIHllYXIiLCBsaW1pdHMgPSBjKDEsNTIpLCAKICAgICAgICAgICAgICAgICAgICAgICAgYnJlYWtzID0gYygxLDI2LDUyKSkgKwogIGxhYnMoeCA9ICJsb2coRmlyc3QgcmF0aW8gb2YgRVdTKSIsIHkgPSBleHByZXNzaW9uKFJbRV0pKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpCmBgYAoKIyMjIENvcnJlbGF0aW9ucyBiZXR3ZWVuICRSX0UkIGFuZCBsb2cgb2YgImZpcnN0IHJhdGlvZWQiIEVXUywgbmVhci16ZXJvcyByZW1vdmVkCgpJdCBsb29rcyBsaWtlIHRoZXJlIG1pZ2h0IGJlIHRyZW5kcyB3aGVuIHRoZSBFV1MgYXJlIGFjdHVhbGx5ICJtb3ZpbmciLCB0aGF0IGlzLCB3aGVuICRcdGV4dHtsb2d9XGxlZnQoXGZyYWN7eCh0KX17eCh0LTEpfSBccmlnaHQpIFxuZSAwJCBvciBhcHByb3hpbWF0ZWx5IHNvLgpTbywgSSdsbCByZW1vdmUgdGhvc2UgdmFsdWVzIHRoYXQgYXJlIGNsb3NlIHRvIHplcm8gYnkgYXNzdW1pbmcgYWxsIHZhbHVlcyBsZXNzIHRoYW4gJHwwLjA1fCQgYXJlIGVzc2VudGlhbGx5IHplcm8uCk5vdyBzb21lIHBhdHRlcm5zIGVtZXJnZSBhcyBleHBlY3RlZC4KCmBgYHtyIGZpZy5oZWlnaHQ9MjAsIGZpZy53aWR0aD04fQpld3NfcmVmZnNfZGlmZmVkICU+JQogIGZpbHRlcihSZWZmIDwgMSkgJT4lCiAgZmlsdGVyKHJvdW5kKGxvZ19yYXRpb192YWx1ZSwxKSAhPSAwKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSBsb2dfcmF0aW9fdmFsdWUsIHkgPSBSZWZmKSkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG9yID0gbHVicmlkYXRlOjp3ZWVrKGRhdGUpKSwgYWxwaGEgPSAwLjcpICsKICBzdGF0X3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBsaW5ldHlwZSA9IDEsIGNvbG9yID0gImJsYWNrIikgKwogIGZhY2V0X3dyYXAoZXdzfnJlZ2lvbiwgc2NhbGVzID0gImZyZWUiLCBuY29sID0gNCkgKwogIHNjYWxlX2NvbG9yX3ZpcmlkaXNfYyhuYW1lID0gIldlZWsgb2YgeWVhciIsIGxpbWl0cyA9IGMoMSw1MiksIAogICAgICAgICAgICAgICAgICAgICAgICBicmVha3MgPSBjKDEsMjYsNTIpKSArCiAgbGFicyh4ID0gImxvZyhGaXJzdCByYXRpbyBvZiBFV1MpIiwgeSA9IGV4cHJlc3Npb24oUltFXSkpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAidG9wIikKYGBgCgojIyBDYWxjdWxhdGUgY29ycmVsYXRpb25zIGFmdGVyIHJlbW92aW5nIH4wIEVXUyBkaWZmcwoKYGBge3IgY2FsYy1jb3JzfQpld3NfcmVmZnNfZGlmZmVkICU+JQogIGZpbHRlcihSZWZmIDwgMSkgJT4lCiAgZmlsdGVyKHJvdW5kKGxvZ19yYXRpb192YWx1ZSwxKSAhPSAwKSAlPiUKICBncm91cF9ieShld3MsIHJlZ2lvbikgJT4lCiAgc3VtbWFyaXNlKGNvcnJfd2l0aF9yZWZmID0gY29yKGxvZ19yYXRpb192YWx1ZSwgUmVmZiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHVzZSA9ICJwYWlyd2lzZS5jb21wbGV0ZS5vYnMiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gInNwZWFybWFuIikpICU+JQogIGdncGxvdChhZXMoeCA9IGV3cywgeSA9IGNvcnJfd2l0aF9yZWZmLCBmaWxsID0gcmVnaW9uKSkgKwogIGdlb21fY29sKHBvc2l0aW9uID0gcG9zaXRpb25fZG9kZ2UoKSwgd2lkdGggPSAwLjUpICsKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gMCksIGNvbCA9ICJncmV5NDAiKSArCiAgbGFicyh5ID0gIlNwZWFybWFuJ3MgcmFuayBjb3JyZWxhdGlvbiIsIHggPSAiRVdTIikgKwogIGdndGhlbWVzOjpzY2FsZV9maWxsX2NvbG9yYmxpbmQoKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGxpbWl0cyA9IGMoLTAuNiwwLjYpKSArCiAgY29vcmRfZmxpcCgpCgp0ZXN0ZGYgPC0gZXdzX3JlZmZzX2RpZmZlZCAlPiUKICBmaWx0ZXIoUmVmZiA8IDEpICU+JQogIGZpbHRlcihyb3VuZChsb2dfcmF0aW9fdmFsdWUsMSkgIT0gMCkgJT4lCiAgZmlsdGVyKGV3cyA9PSAidmFyaWFuY2UiKQpzdW1tYXJ5KGxtKFJlZmZ+bG9nX3JhdGlvX3ZhbHVlLCBkYXRhID0gdGVzdGRmKSkKCmBgYAoKIyBMb29rIGF0IGxhZyBpbiBkeW5hbWljcwpEdWUgdG8gdGhlIHN0b2NoYXN0aWMgc2Vhc29uYWwgZHluYW1pY3MsIHRoZXJlIG1heSBiZSBhIGxhZyBiZXR3ZWVuICRSX0UkIGFuZCBFV1MgdGhhdCBhcmUgY2FsY3VsYXRlZCBmcm9tIHJhdyBjYXNlIGRhdGEuClRvIGxvb2sgYXQgdGhpcywgSSB3aWxsIHBsb3QgJFJfRSQgYW5kIGVhY2ggRVdTIG9uIHRoZSBzYW1lIHBsb3QsIHdpdGggcGVha3MgaW4gZWFjaCB5ZWFyIGhpZ2hsaWdodGVkLgpUaGVzZSBwbG90cyB3aWxsIGNvbWUgZnJvbSBzaW11bGF0ZWQgc2VyaWVzIHdoZXJlIHdlIGtub3cgdGhhdCB0aGUgJFJfRSQgYW5kIHRoZSBkeW5hbWljcyBhcmUgZXhwbGljaXRseSBsaW5rZWQgKHVubGlrZSAkUl9FJCBlc3RpbWF0ZXMgZnJvbSBwYXJ0aWNsZSBmaWx0ZXJpbmcgd2l0aCB0aGUgcmVhbCBkYXRhKS4KCmBgYHtyIGxvYWQtc2ltc30KYWxsX3NpbV9maWxlcyA8LSBsaXN0LmZpbGVzKCIuLi8uLi9zaW11bGF0aW9ucy8iKQpzaW1faWRzIDwtIGdyZXAoImJvb3RzdHJhcCoiLCBhbGxfc2ltX2ZpbGVzKQpzaW1fZmlsZXMgPC0gYWxsX3NpbV9maWxlc1tzaW1faWRzXQoKYWxsX3NpbXMgPC0gdGliYmxlKCkgICMgZW1wdHkgdGliYmxlIGZvciBzdG9yYWdlCmZvcihkb19maWxlIGluIHNpbV9maWxlcyl7CiAgdG1wX3NpbSA8LSByZWFkUkRTKHBhc3RlMCgiLi4vLi4vc2ltdWxhdGlvbnMvIiwgZG9fZmlsZSkpICU+JQogICAgZmlsdGVyKHNpbSA9PSAxKSAlPiUKICAgIHVubmVzdCgpICU+JQogICAgbXV0YXRlKGNpdHkgPSBzdHJfc3ViKGRvX2ZpbGUsIDE2LCAtNSkpCiAgCiAgYWxsX3NpbXMgPC0gYmluZF9yb3dzKGFsbF9zaW1zLCB0bXBfc2ltKQp9CmBgYAoKYGBge3IgY2FsYy1zaW0tZXdzfQpzaW1fZXdzIDwtIHRpYmJsZSgpICAjIGVtcHR5IHN0b3JhZ2UgZGYKZm9yKGRvX2NpdHkgaW4gdW5pcXVlKGFsbF9zaW1zJGNpdHkpKXsKICB0bXBfZGF0YSA8LSBhbGxfc2ltcyAlPiUKICAgIGZpbHRlcihjaXR5ID09IGRvX2NpdHkpCiAgCiAgdG1wX2Nhc2VzIDwtIHRtcF9kYXRhJHJlcG9ydHMKICAKICB0bXBfc3RhdHMgPC0gc3BhZXJvOjpnZXRfc3RhdHMoCiAgICB4ID0gdG1wX2Nhc2VzLAogICAgY2VudGVyX3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBjZW50ZXJfa2VybmVsID0gInVuaWZvcm0iLCAKICAgIGNlbnRlcl9iYW5kd2lkdGggPSAzMCwgCiAgICBzdGF0X3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBzdGF0X2tlcm5lbCA9ICJ1bmlmb3JtIiwgCiAgICBzdGF0X2JhbmR3aWR0aCA9IDMwLCAKICAgIGxhZyA9IDEsIAogICAgYmFja3dhcmRfb25seSA9IFRSVUUKICApJHN0YXRzCiAgCiAgc2ltX2V3cyA8LSBzaW1fZXdzICU+JQogICAgYmluZF9yb3dzKAogICAgICB0bXBfZGF0YSAlPiUgYmluZF9jb2xzKGFzX3RpYmJsZSh0bXBfc3RhdHMpKQogICAgKQp9CmBgYAoKYGBge3IgcGxvdC1keW5hbWljcywgZmlnLmhlaWdodD00LCBmaWcud2lkdGg9MTB9CnNpbV9ld3MgPC0gc2ltX2V3cyAlPiUKICBtdXRhdGUoY3JpdGljYWwgPSBpZmVsc2UoUkVfc2VhcyA8IDEsIEZBTFNFLCBUUlVFKSkKCmZvcihkb19jaXR5IGluIHVuaXF1ZShzaW1fZXdzJGNpdHkpKXsKICB0bXAgPC0gZmlsdGVyKHNpbV9ld3MsIGNpdHkgPT0gZG9fY2l0eSkKICB4IDwtIHRtcCR0aW1lCiAgeSA8LSB0bXAkUkVfc2VhcwogIHkyIDwtIGxvZyh0bXAkdmFyaWFuY2UpCiAgY29sb3IgPC0gaWZlbHNlKHRtcCRjcml0aWNhbCA9PSBUUlVFLCAicmVkIiwgImJsdWUiKQogIAogIHBhcihtYXIgPSBjKDUsNSwyLDUpKQogIHBsb3QoeCwgeSwgdHlwZSA9ICJuIiwgeWxhYiA9IGV4cHJlc3Npb24oaXRhbGljKFIpW0VdKSwgeGxhYiA9ICJEYXRlIiwgbWFpbiA9IGRvX2NpdHkpCiAgc2VnbWVudHMoeDAgPSB4LCB5MCA9IHksIHgxID0geCtkaWZmKHgpLCB5MSA9IHkrZGlmZih5KSwgY29sID0gY29sb3IpCiAgcGFyKG5ldyA9IFQpCiAgcGxvdCh4LCB5MiwgdHlwZSA9ICJsIiwgYXhlcz1GLCB4bGFiPU5BLCB5bGFiPU5BLCBsd2QgPSAxLjUpCiAgYXhpcyhzaWRlID0gNCkKICBtdGV4dChzaWRlID0gNCwgbGluZSA9IDMsICJMb2cgdmFyaWFuY2UiKQp9Cgpmb3IoZG9fY2l0eSBpbiB1bmlxdWUoc2ltX2V3cyRjaXR5KSl7CiAgdG1wIDwtIGZpbHRlcihzaW1fZXdzLCBjaXR5ID09IGRvX2NpdHkpCiAgeCA8LSB0bXAkdGltZQogIHkgPC0gdG1wJFJFX3NlYXMKICB5MiA8LSB0bXAkYXV0b2NvcnJlbGF0aW9uCiAgY29sb3IgPC0gaWZlbHNlKHRtcCRjcml0aWNhbCA9PSBUUlVFLCAicmVkIiwgImJsdWUiKQogIAogIHBhcihtYXIgPSBjKDUsNSwyLDUpKQogIHBsb3QoeCwgeSwgdHlwZSA9ICJuIiwgeWxhYiA9IGV4cHJlc3Npb24oaXRhbGljKFIpW0VdKSwgeGxhYiA9ICJEYXRlIiwgbWFpbiA9IGRvX2NpdHkpCiAgc2VnbWVudHMoeDAgPSB4LCB5MCA9IHksIHgxID0geCtkaWZmKHgpLCB5MSA9IHkrZGlmZih5KSwgY29sID0gY29sb3IpCiAgcGFyKG5ldyA9IFQpCiAgcGxvdCh4LCB5MiwgdHlwZSA9ICJsIiwgYXhlcz1GLCB4bGFiPU5BLCB5bGFiPU5BLCBsd2QgPSAxLjUpCiAgYXhpcyhzaWRlID0gNCkKICBtdGV4dChzaWRlID0gNCwgbGluZSA9IDMsICJBdXRvY29ycmVsYXRpb24iKQp9Cgpmb3IoZG9fY2l0eSBpbiB1bmlxdWUoc2ltX2V3cyRjaXR5KSl7CiAgdG1wIDwtIGZpbHRlcihzaW1fZXdzLCBjaXR5ID09IGRvX2NpdHkpCiAgeCA8LSB0bXAkdGltZQogIHkgPC0gdG1wJFJFX3NlYXMKICB5MiA8LSBsb2codG1wJGRlY2F5X3RpbWUgKyAxKQogIGNvbG9yIDwtIGlmZWxzZSh0bXAkY3JpdGljYWwgPT0gVFJVRSwgInJlZCIsICJibHVlIikKICAKICBwYXIobWFyID0gYyg1LDUsMiw1KSkKICBwbG90KHgsIHksIHR5cGUgPSAibiIsIHlsYWIgPSBleHByZXNzaW9uKGl0YWxpYyhSKVtFXSksIHhsYWIgPSAiRGF0ZSIsIG1haW4gPSBkb19jaXR5KQogIHNlZ21lbnRzKHgwID0geCwgeTAgPSB5LCB4MSA9IHgrZGlmZih4KSwgeTEgPSB5K2RpZmYoeSksIGNvbCA9IGNvbG9yKQogIHBhcihuZXcgPSBUKQogIHBsb3QoeCwgeTIsIHR5cGUgPSAibCIsIGF4ZXM9RiwgeGxhYj1OQSwgeWxhYj1OQSwgbHdkID0gMS41KQogIGF4aXMoc2lkZSA9IDQpCiAgbXRleHQoc2lkZSA9IDQsIGxpbmUgPSAzLCAiTG9nIGRlY2F5IHRpbWUgKCsxKSIpCn0KCgpgYGA=